#### prepare a function ####
issue.var.name <- c("a.article9", "b.defense", "c.revisionism", "d.women", 
                    "e.gay", "f.immigrant", "g.growth", "h.tax")

# function to conduct F-tests for interactions between attributes
interaction.F.test <- function(variable.1, variable.2, data, choice) {
  if (choice == TRUE) {
    lm.model.1 <- lm(choice ~ a.article9 + b.defense + c.revisionism + d.women + 
                       e.gay + f.immigrant + g.growth + h.tax, data = data)
  } else {
    lm.model.1 <- lm(rating ~ a.article9 + b.defense + c.revisionism + d.women + 
                       e.gay + f.immigrant + g.growth + h.tax, data = data)
  }
  if (choice == TRUE) {
    lm.formula <- as.formula(paste0("choice ~ a.article9 + b.defense + c.revisionism + d.women + 
                                      e.gay + f.immigrant + g.growth + h.tax + ", 
                                    variable.1, " : ", variable.2))
  } else {
    lm.formula <- as.formula(paste0("rating ~ a.article9 + b.defense + c.revisionism + d.women + 
                                      e.gay + f.immigrant + g.growth + h.tax + ", 
                                    variable.1, " : ", variable.2))
  }
  lm.model.2 <- lm(lm.formula, data = data)
  anova(lm.model.1, lm.model.2)$"Pr(>F)"[2]
}

#### load data ####
task.data <- read.csv("task_data.csv")

# reorder the levels of attribute variables
position.label <- c("Agree", "Neither", "Disagree")
task.data$a.article9 <- factor(task.data$a.article9, levels = position.label)
task.data$b.defense <- factor(task.data$b.defense, levels = position.label)
task.data$c.revisionism <- factor(task.data$c.revisionism, levels = position.label)
task.data$d.women <- factor(task.data$d.women, levels = position.label)
task.data$e.gay <- factor(task.data$e.gay, levels = position.label)
task.data$f.immigrant <- factor(task.data$f.immigrant, levels = position.label)
task.data$g.growth <- factor(task.data$g.growth, levels = position.label)
task.data$h.tax <- factor(task.data$h.tax, levels = position.label)

#### analysis ####
## conduct F-tests
interaction.F.test.p.choice <- 
  interaction.F.test.p.rating <- matrix(NA, 8, 8)
for (i in 1:8) {
  for (j in 1:8) {
    if (i == j) break
    interaction.F.test.p.choice[i, j] <- 
      interaction.F.test(issue.var.name[i], issue.var.name[j], 
                         data = task.data, choice = TRUE)
    interaction.F.test.p.rating[i, j] <- 
      interaction.F.test(issue.var.name[i], issue.var.name[j], 
                         data = task.data, choice = FALSE)
  }
}

# Table A.12
round(interaction.F.test.p.choice, 3)
round(interaction.F.test.p.rating, 3)

# number of pairs passing the F-test
sum(interaction.F.test.p.choice < 0.05, na.rm = TRUE)
sum(interaction.F.test.p.rating < 0.05, na.rm = TRUE)

# p-value adjustment by the Benjamini-Hochberg method
round(matrix(p.adjust(interaction.F.test.p.choice, method = "BH"), 8, 8), 3)
round(matrix(p.adjust(interaction.F.test.p.rating, method = "BH"), 8, 8), 3)